#ifndef _CALC_H_
#define _CALC_H_
//#ifndef _GNU_SOURCE
//#define _GNU_SOURCE
//#include<iostream>
#include <sched.h>
#include <pthread.h>
//using namespace std;
#ifdef __cplusplus
extern "C"{	
#endif

typedef struct {
	int* 			rowArray;
	const int* 		rowOffset;
	int 			rowArraySize;
	const int* 		columnIndice;
	const double* 	S;
	const double*	valueNormalMatrix;
	double* 		Id;
}TaskMatrixInfoA;

typedef struct {
	const double * 	valueSpiceMatrix;//N*2
	const int* 		rowOffset; 		//n+1
	const int* 		columnIndice; 	//N

	double* 		A;	//N
	double*		 	S;	//nx1
	double*			R;	//nx1
	double* 		H;	//nx1
	const double* 	D;	//N*2
	double* 		IC;	//nx1
	double* 		IG;	//nx1
	double 		alpha;
	int * 			rowArray; //kx1
	int 			rowArraySize;
	void* 			hdl;  //? what's this?
}TaskMatrixInfoB;


typedef struct {
	int n1;
	int n2;
	TaskMatrixInfoA ** listDataList;
}MY_ARGSA;
#define THREAD_NUM 32

void* mysuma(void* args){
	MY_ARGSA* my_args = (MY_ARGSA*) args;
	int n1 = my_args->n1;
	int n2 = my_args->n2;
	TaskMatrixInfoA **listDataList = my_args->listDataList;

	for (int i = n1 ;i < n2; i=i+THREAD_NUM)
	{	
		int row_array_size= listDataList[i]->rowArraySize;	
		TaskMatrixInfoA * list_Data_List=listDataList[i];
		const double *valueNormalMatrix = list_Data_List->valueNormalMatrix;
		const double *S = list_Data_List->S;
		const int *columnIndice = list_Data_List->columnIndice;
		for (unsigned int it = 0; it < row_array_size; ++it	)
		{
			const int node = list_Data_List->rowArray[it];
			int jstart = list_Data_List->rowOffset[node];
			int jend = list_Data_List->rowOffset[node + 1];
			double tmp=0;
			int nnz = jend - jstart;			
			if (nnz <= 4) {
				switch (nnz) {
					case 1:
						tmp += valueNormalMatrix[jstart]*S[columnIndice[jstart]];
						break;
					case 2:
						tmp += valueNormalMatrix[jstart]*S[columnIndice[jstart]]
								+valueNormalMatrix[jstart+1]*S[columnIndice[jstart+1]];
						break;
					case 3:				
						tmp += valueNormalMatrix[jstart]*S[columnIndice[jstart]]
							+valueNormalMatrix[jstart+1]*S[columnIndice[jstart+1]]
							+valueNormalMatrix[jstart+2]*S[columnIndice[jstart+2]];
						break;
					case 4:
						tmp += valueNormalMatrix[jstart]*S[columnIndice[jstart]]
							+valueNormalMatrix[jstart+1]*S[columnIndice[jstart+1]]
							+valueNormalMatrix[jstart+2]*S[columnIndice[jstart+2]]
							+valueNormalMatrix[jstart+3]*S[columnIndice[jstart+3]];
				}
			} else {
				int m = nnz%4;
				switch (m) {
					case 1:
						tmp += valueNormalMatrix[jstart]*S[columnIndice[jstart]];
						break;
					case 2:
						tmp += valueNormalMatrix[jstart]*S[columnIndice[jstart]]
								+valueNormalMatrix[jstart+1]*S[columnIndice[jstart+1]];
						break;
					case 3:				
						tmp += valueNormalMatrix[jstart]*S[columnIndice[jstart]]
							+valueNormalMatrix[jstart+1]*S[columnIndice[jstart+1]]
							+valueNormalMatrix[jstart+2]*S[columnIndice[jstart+2]];
						break;
				}
				for (int j = jstart+m;j < jend;j+=4)
				{
					tmp += valueNormalMatrix[j]*S[columnIndice[j]]
						 +valueNormalMatrix[j+1]*S[columnIndice[j+1]]
						 +valueNormalMatrix[j+2]*S[columnIndice[j+2]]
						 +valueNormalMatrix[j+3]*S[columnIndice[j+3]];
				}
			}
			list_Data_List->Id[node]+=tmp;
		}
	}
	return NULL;
}

typedef struct {
	int n1;
	int n2;
	TaskMatrixInfoB ** listDataList;
}MY_ARGSB;


void* mysumb(void* args){
	MY_ARGSB* my_args = (MY_ARGSB*) args;
	int n1 = my_args->n1;
	int n2 = my_args->n2;
	TaskMatrixInfoB **listDataList = my_args->listDataList;

	for (int i = n1;i<n2;i=i+THREAD_NUM){
		TaskMatrixInfoB * list_Data_List = listDataList[i];
		int row_array_size = list_Data_List->rowArraySize;
		const double alpha = list_Data_List->alpha;
		double* A = list_Data_List->A;
		const double* S = list_Data_List->S;
		const int* cI=list_Data_List->columnIndice;
		const double* vSM=list_Data_List->valueSpiceMatrix;

		for(unsigned int it = 0;it < row_array_size;++it)
		{
			int node = list_Data_List->rowArray[it];
			const int k1 = node <<1;
									
			unsigned int  jstart = list_Data_List->rowOffset[node];
			unsigned int  jend = list_Data_List->rowOffset[node + 1];
			unsigned nnz = jend-jstart; 
			unsigned int k = jstart<<1;

			double ig_tmp =0;
			double ic_tmp =0;
			if(nnz<=4){
				switch(nnz){
					case 1:{
						ig_tmp += vSM[k] *S[cI[jstart]];
						ic_tmp += vSM[k+1]*S[cI[jstart]];
						A[jstart]=vSM[k]+vSM[k+1]*alpha;
						break;
					}
					case 2:{
						ig_tmp += vSM[k]  *S[cI[jstart]]+vSM[k+2]*S[cI[jstart+1]];
						ic_tmp += vSM[k+1]*S[cI[jstart]]+vSM[k+3]*S[cI[jstart+1]];
						A[jstart]  =vSM[k]  +vSM[k+1]*alpha;
						A[jstart+1]=vSM[k+2]+vSM[k+3]*alpha;
						break;
					}
					case 3:{
						ig_tmp += vSM[k]  *S[cI[jstart]]+vSM[k+2]*S[cI[jstart+1]]+vSM[k+4]*S[cI[jstart+2]];
						ic_tmp += vSM[k+1]*S[cI[jstart]]+vSM[k+3]*S[cI[jstart+1]]+vSM[k+5]*S[cI[jstart+2]];
						A[jstart]  =vSM[k]  +vSM[k+1]*alpha;
						A[jstart+1]=vSM[k+2]+vSM[k+3]*alpha;
						A[jstart+2]=vSM[k+4]+vSM[k+5]*alpha;
						break;
					}
					case 4:{
						ig_tmp += vSM[k]  *S[cI[jstart]]+vSM[k+2]*S[cI[jstart+1]]+vSM[k+4]*S[cI[jstart+2]]+vSM[k+6]*S[cI[jstart+3]];
						ic_tmp += vSM[k+1]*S[cI[jstart]]+vSM[k+3]*S[cI[jstart+1]]+vSM[k+5]*S[cI[jstart+2]]+vSM[k+7]*S[cI[jstart+3]];
						A[jstart]  =vSM[k]  +vSM[k+1]*alpha;
						A[jstart+1]=vSM[k+2]+vSM[k+3]*alpha;
						A[jstart+2]=vSM[k+4]+vSM[k+5]*alpha;
						A[jstart+3]=vSM[k+6]+vSM[k+7]*alpha;
					}
				}
				
			}else{
				int m = nnz%4;
				switch(m){
					case 1:{
						ig_tmp += vSM[k] *S[cI[jstart]];
						ic_tmp += vSM[k+1]*S[cI[jstart]];
						A[jstart]=vSM[k]+vSM[k+1]*alpha;
						break;
					}
					case 2:{
						ig_tmp += vSM[k]  *S[cI[jstart]]+vSM[k+2]*S[cI[jstart+1]];
						ic_tmp += vSM[k+1]*S[cI[jstart]]+vSM[k+3]*S[cI[jstart+1]];
						A[jstart]  =vSM[k]  +vSM[k+1]*alpha;
						A[jstart+1]=vSM[k+2]+vSM[k+3]*alpha;
						break;
					}
					case 3:{
						ig_tmp += vSM[k]  *S[cI[jstart]]+vSM[k+2]*S[cI[jstart+1]]+vSM[k+4]*S[cI[jstart+2]];
						ic_tmp += vSM[k+1]*S[cI[jstart]]+vSM[k+3]*S[cI[jstart+1]]+vSM[k+5]*S[cI[jstart+2]];
						A[jstart]  =vSM[k]  +vSM[k+1]*alpha;
						A[jstart+1]=vSM[k+2]+vSM[k+3]*alpha;
						A[jstart+2]=vSM[k+4]+vSM[k+5]*alpha;
						break;
					}
				}
				for (int j = jstart+m;j<jend;j+=4){
					int k=j<<1;
					ig_tmp += vSM[k]  *S[cI[j]]+vSM[k+2]*S[cI[j+1]]+vSM[k+4]*S[cI[j+2]]+vSM[k+6]*S[cI[j+3]];
					ic_tmp += vSM[k+1]*S[cI[j]]+vSM[k+3]*S[cI[j+1]]+vSM[k+5]*S[cI[j+2]]+vSM[k+7]*S[cI[j+3]];
					A[j]  =vSM[k]  +vSM[k+1]*alpha;
					A[j+1]=vSM[k+2]+vSM[k+3]*alpha;
					A[j+2]=vSM[k+4]+vSM[k+5]*alpha;
					A[j+3]=vSM[k+6]+vSM[k+7]*alpha;
				}
			}

			list_Data_List->IG[node] += ig_tmp;
 			list_Data_List->IC[node] += ic_tmp;
			list_Data_List->R[node] = list_Data_List->D[k1] -ig_tmp;//current; 
			list_Data_List->H[node] = list_Data_List->D[k1+1] - ic_tmp;//charge;
		} 
	}
	return NULL;
}

void matrix_calc_taskA(TaskMatrixInfoA** listDataList,int N)
{
	int thread_num = THREAD_NUM;
	pthread_t th[thread_num];

	int i;
	MY_ARGSA args[thread_num] ;
	for(i =0;i<thread_num;i++){
		args[i].n1 = i;
		args[i].n2 = N;
		args[i].listDataList = listDataList;
	}

	//设置亲和性
	cpu_set_t mask[20];
	pthread_attr_t attr[20];
	for(i =0;i<20;i++){		
		CPU_ZERO(&mask[i]);
		CPU_SET(i,&mask[i]);
		CPU_SET(i+20,&mask[i]);
		pthread_attr_init(&attr[i]);
		pthread_attr_setaffinity_np(&attr[i], sizeof(cpu_set_t), &mask[i]);
	}
	sched_setaffinity(0, sizeof(cpu_set_t), &mask[19]);
	//产生线程	
	for(i =0;i<thread_num;i++){
		pthread_create(&th[i],&attr[i/2],mysuma,&args[i]);		
	}

	for(i =0;i<thread_num;i++){
		pthread_join(th[i],NULL);
	}
}

//Task B interface
void matrix_calc_taskB(TaskMatrixInfoB** listDataList,int N)
{
	int thread_num = THREAD_NUM;
	pthread_t th[thread_num];

	int i;
	MY_ARGSB args[thread_num] ;
	for(i =0;i<thread_num;i++){
		args[i].n1 = i;
		args[i].n2 = N;
		args[i].listDataList = listDataList;
	}

	cpu_set_t mask[20];
	pthread_attr_t attr[20];
	for(i =0;i<20;i++){		
		CPU_ZERO(&mask[i]);
		CPU_SET(i,&mask[i]);
		CPU_SET(i+20,&mask[i]);
		pthread_attr_init(&attr[i]);
		pthread_attr_setaffinity_np(&attr[i], sizeof(cpu_set_t), &mask[i]);
	}
	sched_setaffinity(0, sizeof(cpu_set_t), &mask[19]);
	//产生线程	
	for(i =0;i<thread_num;i++){
		pthread_create(&th[i],&attr[i/2],mysumb,&args[i]);		
	}

	for(i =0;i<thread_num;i++){
		pthread_join(th[i],NULL);
	}
	
}

#ifdef __cplusplus
}
#endif

#endif
